Load Data for Point Base Map: Mobile Home Data


source: https://hifld-geoplatform.opendata.arcgis.com/datasets/mobile-home-parks

Using data from the US Department of Homeland Security, I have plotted locations of trailer parks across New York State. It appears most trailer parks have fewer than 50 trailers and are unsurprisingly concentrated Upstate.


# Marker Point Map

# Load mobile home data
pal <- colorFactor(c("mediumturquoise",  "purple" , "cadetblue4" ), domain = c("LARGE (>100)", "MEDIUM (51-100)","SMALL (<50)"))

basemap_points <-leaflet(mobile_home_ny) %>%
  # Basemaps
  addTiles()%>%
  addProviderTiles("Esri.WorldPhysical", group = "WorldPhysical") %>%
  addCircleMarkers(group = mobile_home_ny$SIZE, stroke = F, color = ~pal(SIZE), weight = 50,
                   fillOpacity = 0.5) %>%
    addLayersControl(
    baseGroups = c("OSM (default)", "WorldPhysical"),
    overlayGroups = c(levels(mobile_home_ny$SIZE)),
    options = layersControlOptions(collapsed = FALSE))%>%
   addLegend("bottomright", pal = pal, values = ~SIZE,
    title = "Size of Mobile Home Community")
## Assuming "LONGITUDE" and "LATITUDE" are longitude and latitude, respectively
basemap_points %>% 
  fitBounds(-78, 45, -72, 40) # Fit to NY


# Data from : https://hifld-geoplatform.opendata.arcgis.com/datasets/public-transit-routes/data?geometry=-104.352%2C31.109%2C-57.023%2C43.309

Public Transit Data from DHS


This plot shows the public transit systems and the system authority under which the transit system runs under. Again, most of these transit systems are concentrated in New York City, with a very small transit line in Buffalo.

# LINE PLOT

# Load Shape File 
public_transit <- readOGR("Public_Transit_Routes.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "G:\My Drive\Fall 2019\Mini 1\RShiny for Operations Management\HW\HW_2\HW2_nmurray\Public_Transit_Routes.shp", layer: "Public_Transit_Routes"
## with 316 features
## It has 17 fields
## Integer64 fields read as strings:  OBJECTID ID UACODE
# Limit to NY Transit Systems
ny_lines <- c("METRO-NORTH COMMUTER RAILROAD COMPANY, DBA: MTA METRO-NORTH RAILROAD",  
              "MTA LONG ISLAND RAIL ROAD",                                                                  
              "MTA NEW YORK CITY TRANSIT",                                                                  
              "NIAGARA FRONTIER TRANSPORTATION AUTHORITY",                                                  
              "PORT AUTHORITY OF NEW YORK AND NEW JERSEY",                                                  
              "PORT AUTHORITY TRANS-HUDSON CORPORATION",                                                    
              "STATEN ISLAND RAPID TRANSIT OPERATING AUTHORITY, DBA: MTA STATEN ISLAND RAILWAY")

public_transit_ne <- subset(public_transit, public_transit$SYS_AGENCY %in% ny_lines)

#Create Color Pallette

cols <- colorRampPalette(brewer.pal(8,"Dark2"))(length(ny_lines))
pal2 <-colorFactor(cols, domain = ny_lines)

leaflet(data = public_transit_ne) %>%
  # Basemaps
  addTiles() %>%
  addProviderTiles("Esri.WorldGrayCanvas", group = "WorldGrayCanvas") %>%
  addProviderTiles("Esri.WorldTerrain", group = "WorldTerrain") %>%
  addPolylines(group = public_transit_ne$SYS_AGENCY, color = ~pal2(SYS_AGENCY))  %>%
  # Layers control
  addLayersControl(
    baseGroups = c("WorldGrayCanvas", "WorldTerrain"),
    # overlayGroups = ne_lines,
    # position = "bottomleft",
    options = layersControlOptions(collapsed = FALSE)) %>%
   addLegend("topright", pal = pal2, values = ~SYS_AGENCY,
    title = " New York Public Transit Lines By System")  %>%
    fitBounds(-79, 48, -72, 39) # Fit to Northeast

Polygon Map with 2 Layers


Given that we know there are more mobiles homes Upstate and more transit lines downstate, one might infer that Downstate is where the population center is.


If Downstate has more transit lines, do they use fewer cars or does the difference in population outweight this? Based on this map, we conclude that Downstate has a much larger population and the different of population means more cares, irrespective of public transit.

Prep the data and add shapefile


Data source: https://www.census.gov/data.html

# Source: https://walkerke.github.io/tidycensus/reference/load_variables.html
# Get the data from the census

v17 <- load_variables(2017, "acs5", cache = TRUE)

# Vehicle Ownership 
ny_data_car <- get_acs(geography = "county", 
                       variables = c("B08541_006"),
                       state = 36,
                       geometry = FALSE)
## Getting data from the 2012-2016 5-year ACS
# Population
ny_data_pop <- get_acs(geography = "county", 
                       variables = c("B01003_001"),
                       state = 36,
                       geometry = FALSE)
## Getting data from the 2012-2016 5-year ACS
# Get NYS county shapefile 
ny_counties <- counties(state = 36)


# Merge Car and Population data files in a single file 

poly_map_data <- merge(ny_data_car, ny_data_pop, by = c("GEOID", "NAME"), sort = FALSE)

ny_counties@data<- merge(ny_counties@data, poly_map_data, by = c("GEOID"), sort = FALSE)

#Merge Car and Shapefile data
cars_ny <- merge(ny_counties, ny_data_car, by = c("GEOID"))


# Recall Varible x is car data and variable.y is population data 

Polygon Map with Number of Cars and Number of People

# Palette for Number of Cars 
pal_car <- colorNumeric(
  palette = "Reds",
  na.color = "#808080",
  domain = ny_counties$estimate.x)  # of Cars 

# Palette for Population
pal_ppl <- colorNumeric(
  palette = "Purples",
  na.color = "#808080",
  domain = ny_counties$estimate.y)  # population 

basemap_poly <- leaflet(ny_counties) %>%
  addTiles() %>%
  addProviderTiles("Esri.WorldGrayCanvas", group = "WorldGrayCanvas") %>%
  addProviderTiles("CartoDB.Positron", group = "Positron") %>%
    fitBounds(-79, 45, -72, 39) # Fit to Northeast

basemap_poly %>%
  addPolygons(
        color = ~pal_car(estimate.x),
        stroke = FALSE,
        weight = 2,
        opacity = 1,
        dashArray = "3",
        fillOpacity = 0.7, 
        group = "Number of Cars") %>%
            addLegend(group = "Number of Cars", "bottomright", 
                      pal = pal_car, values = ny_counties$estimate.x, 
                      data = getMapData(basemap_poly),title = "Number of Cars") %>%
  addPolygons(
        color = ~pal_ppl(estimate.y),
        stroke = FALSE,
        weight = 2,
        opacity = 1,
        dashArray = "3",
        fillOpacity = 1, 
        group = "Number of People")%>%
            addLegend(group = "Number of People", "bottomright", 
                      pal = pal_ppl, values = ny_counties$estimate.y, 
                      data = getMapData(basemap_poly), title = "Number of People") %>%
  addLayersControl(
      baseGroups = c("WorldGrayCanvas", "Positron"),
      overlayGroups = c("Number of Cars", "Number of People"),
      options = layersControlOptions(collapsed = FALSE)) %>% 
  hideGroup("Number of People")